^DEPARTMENT  OF  DEFEN cT 

DEFENCE  SCIENCE  &  TECHNOLOGY  ORGANISATION 


DSTO 


Determining  Beam  Bending 
Distribution  Using  Dynamic 
Information 

Frank  G.  Polanco 
DSTO-RR-0226 


DISTRIBUTION  STATEMENT  A: 
Approved  for  Public  Release  - 
Distribution  Unlimited 


20020226  084 


Determining  Beam  Bending  Distribution  Using 
Dynamic  Information 

Frank  G.  Polanco 


Airframes  and  Engines  Division 
Aeronautical  and  Maritime  Research  Laboratory 


DSTO-RR-0226 


ABSTRACT 


As  a  first  approxinration,  a  helicopter  rotor  blade  may  be  modelled  as  a 
cantilever  beam.  Given  the  initial  deformation  of  this  beam,  and  using  either 
strain  or  acceleration  at  one  location  along  the  beam,  we  can  determine  the 
load  distribution  along  the  entire  beam.  We  consider  load  distributions  that 
can  vary  spatially,  but  are  constant  in  time  (except  for  the  initial  step  input). 
In  the  solution  we  neglect  the  effects  of  both  aerodynamic  and  mechanical 
damping.  The  separation  of  variables  technique  leads  to  a  solution  in  terms  of 
the  beam’s  natural  modes.  The  loading  distribution  is  decomposed  in  terms 
of  these  modes.  A  finite  element  simulation  of  the  beam’s  response  to  a  cubic 
load  distribution  verifies  that  this  load  prediction  is  possible.  We  demonstrate 
that  the  higher  modes  of  the  load  prediction  are  unstable  when  noise  is  present 
in  the  measurements,  but  that  the  lower  modes  are  robust.  If  the  initial  beam 
deformation  is  unknown,  then  additional  (strain  or  vibration)  measurement 
locations  may  be  substituted  for  the  unknown  initial  deformation. 
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Determining  Beam  Bending  Distribution  Using 
Dynamic  Information 


EXECUTIVE  SUMMARY 

In  service,  helicopter  rotor  blades  undergo  a  load  distribution  that  is  currently  not 
measured;  yet  this  loading  distribution  determines  the  fatigue  life  of  a  blade.  Given 
that  the  rotor  blades  are  critical  components  in  helicopters,  this  lack  of  information  has 
safety  implications.  Currently,  to  maintain  a  conservative  safety  factor  rotor  blades  are 
necessarily  retired  well  before  their  usage  life  has  been  consumed.  Accurately  predicting 
the  load  distribution  on  a  helicopter  rotor  blade  has  the  potential  to  both  increase  safety 
and  reduce  operating  costs.  In  this  report,  we  use  the  vibration  of  the  rotor  blade  to 
determine  the  loading  distribution  along  the  blade. 

The  rotor  blade  is  modelled  as  a  simple  cantilever  beam  (although  more  general  end 
constraints  are  possible),  and  both  aerodynamic  and  mechanical  damping  are  ignored.  We 
begin  our  load  determination  analysis  with  the  beam  already  vibrating  due  to  a  previous 
excitation.  Hence,  when  we  begin  our  analysis  the  beam  already  has  both  an  initial 
deformation  and  velocity.  The  unknown  load  distribution  is  then  applied  as  a  step  function, 
and  remains  temporally  constant  after  its  application.  In  order  to  determine  the  loading 
distribution,  we  require  the  beam’s  vibration  response  to  this  step  load.  Either  strains 
or  accelerations  may  be  used  as  the  beam’s  vibration  response.  We  show  that  given  the 
beam’s  initial  deformation  and  vibration  response  at  only  one  location,  it  is  possible  to 
predict  the  load  distribution  along  the  entire  blade. 

We  use  the  separation  of  variables  method  to  write  the  solution  in  terms  of  the  blade’s 
natural  modes.  We  then  solve  for  the  load  distribution  in  terms  of  these  natural  modes. 
It  is  not  possible  to  determine  both  the  initial  deformation  and  loading  distribution  given 
only  the  vibration  response  at  one  location. 

A  finite  element  (FE)  model  of  the  beam’s  response  to  a  step  load  is  used  to  simulate 
the  strain  near  the  fixed  end  of  the  cantilever  beam.  We  use  these  simulated  strain  data  to 
predict  the  load  distribution  along  the  FE  beam.  The  simulations  verify  the  load  prediction 
capability  outlined,  but  demonstrate  the  prediction’s  susceptibility  to  measurement  noise 
(especially  in  predicting  higher  frequency  modes). 

If  the  initial  beam  deformation  is  unknown,  then  additional  (strain  or  vibration)  mea¬ 
surement  locations  may  be  substituted  for  knowledge  of  the  initial  deformation. 
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Notation 

Roman  Symbols  Definition 

a  Fourier  cosine  coefficient  (p.  14) 

A  constant  from  fourth  order  differential  equation  (p.  5) 

Aj  constant  from  jth  second  order  differential  equation  (p.  10) 

b  Fourier  sine  coefficient  (p.  14) 

B  constant  from  fourth  order  differential  equation  (p.  5) 

Bj  constant  from  jth  second  order  differential  equation  (p.  10) 

C  constant  from  fourth  order  differential  equation  (p.  5) 

D  constant  from  fourth  order  differential  equation  (p.  5) 

At  time  step  between  temporal  samples  (p.  18) 

E  Young’s  modulus  (p.  3) 

F  axial  load  on  a  beam  (p.  27) 

/  distributed  force  acting  along  beam  (p.  3) 

h  hat  function  (sum  of  two  step  functions)  (p.  12) 

i  dummy  variable  (sometimes  used  as  count  for  discrete  load)  (p.  9) 

i  dummy  variable  (sometimes  used  as  count  for  Fourier  series)  (p.  14) 

j  dummy  variable  (sometimes  used  as  count  for  series  expansion)  (p.  8) 

I  beam’s  cross-sectional  area  moment  of  inertia  (p.  3) 

K  number  of  load  segments  (p.  9) 

I  beam  length  (p.  3) 

m  beam  mass  (p.  3) 

n  number  of  mode  shapes  in  approximating  solution  (p.  15) 

P  loading  period  (p.  13) 

t  temporal  variable  (p.  3) 

ti  counting  variable  (for  the  evaluation  of  relative  errors)  (p.  18) 

tfin  final  count  (number  of  time  evaluations)  (p.  18) 

T  temporal  variable  in  separation  of  variables  technique  (p.  5) 

V  transverse  velocity  of  beam  (p.  4) 

X  horizontal  spatial  variable  (distance  along  beam)  (p.  3) 

y  transverse  deformation  of  beam  (p.  3) 

Y  transverse  deformation  variable  in  separation  of  variables  (p.  5) 


notation  continued  on  next  page  . . . 
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. . .  notation  continued  from  previous  page 


Greek  Symbols 

Definition 

a 

hyperbolic  constant:  function  of  /3 

(p.  5) 

/? 

separation  constant:  from  the  separation  of  variables 

(p.  5) 

7 

represents  series  from  beam  acceleration 

(p.  11) 

6 

delta  (or  Kronecker  delta)  function 

(p.  9) 

4> 

coefficient  for  series  expansion  of  discrete  load 

(p.  8) 

$ 

magnitude  of  discrete  load 

(p.  9) 

V; 

coefficient  of  series  expansion  of  vertical  deformation 

(p.  8) 

a 

initial  beam  deformation  decomposed  into  modes 

(p.  14) 

r 

initial  beam  velocity  decomposed  into  modes 

(p.  14) 

location  of  discrete  load 

(p.  9) 

accelerometer  or  strain  gauge  location  along  beam 

(p.  11) 

U) 

loading  frequency 

(p.  14) 

Miscellaneous  Symbols 

(in  the  following  notation  a  is  a  dummy  variable) 


Oj 

denotes  the  jth  solution  of  the  separation  of  variables 

(P- 

5) 

Oo 

initial  condition 

(P- 

4) 

a, 

dimensional  variable 

(P- 

3) 

a' 

derivative  of  o,  with  respect  to  space,  if  a  = 

-  a(x)  then  a'  =  da/dx 

(P- 

5) 

a, 

derivative  of  a,  with  respect  to  time,  if  o  = 

a{t)  then  a,  =  da/dt 

(P- 

5) 

nth  derivative  of  a,  if  a  =  a{x)  then  = 

dda/dxt^ 

(P- 

5) 

0(0") 

of  order  a" 

(P- 

7) 

a 

approximation  of  a 

(P- 

7) 
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1  Introduction 

Helicopter  rotor  blades  are  considered  critical  components;  yet  the  loading  these  blades 
undergo  in  actual  service  is  largely  unknown.  Determining  the  loading  distribution  history 
of  a  helicopter  rotor  blade  has  the  potential  to  lead  to  increased  safety.  Additionally, 
because  rotor  blades  are  currently  retired  at  a  specified  component  retirement  time,  these 
blades  may  very  well  have  a  significant  amount  of  life  (at  an  acceptable  level  of  reliability) 
remaining.  The  accurate  determination  of  the  in-flight  loading  distribution  on  a  rotor 
blade  therefore  has  the  potential  to  save  the  helicopter  operator  unnecessary  replacement 
costs. 

The  approach  we  use  is  analogous  to  determining  the  load  around  the  rim  of  a  glass 
given  the  vibration  response  of  the  glass.  The  glass  will  ring  differently  depending  on  the 
loading  distribution  around  the  rim.  For  a  rotor  blade  this  “ringing”  is  unique. 

In  this  report  we  develop  a  method  to  determine  the  load  distribution  on  a  cantilever 
beam  (a  first  order  approximation  of  a  helicopter  rotor  blade).  We  determine  this  load 
distribution  using  either  strain  or  acceleration  measurements  at  only  one  location  along 
the  beam. 

Unfortunately  we  can  only  determine  this  loading  distribution  if  we’re  given  the  beam’s 
initial  deformation.  Although  we  model  a  rotor  blade  only  as  a  cantilever  beam,  a  similar 
analysis  may  be  applied  to  more  generalised  end  conditions  (for  example,  an  end  constraint 
that  is  a  combination  of  a  spring  and  pin  joint). 

In  the  remainder  of  this  introduction  (§  1.1)  we  will  review  the  literature  for  similar 
problems.  In  §  2  we  develop  the  partial  differential  equation  (PDE)  governing  the  vibration 
of  a  cantilever  beam.  In  §  3  we  solve  this  PDE  using  the  well  known  separation  of  variables 
technique.  In  §  4  we  develop  the  vibration  solution  for  different  loading  distributions,  which 
include  impulse  loads,  step  loads,  and  a  distributed  load  in  terms  of  the  beam’s  natural 
modes.  In  §  5  we  perform  a  simulation,  which  verifies  the  theoretical  solutions  we  develop. 
We  then  discuss  some  of  the  assumptions  and  findings  in  §  6  before  concluding  the  report 
in  §  7. 


1.1  Review  of  Beam  Vibration 

Barcilon  [3]  determined  the  elastic  properties  (that  is,  stiffness  and  mass  distribution) 
of  a  vibrating  body  using  measurement  data.  As  an  example  a  discretised  beam  was 
analysed,  one  end  was  free  and  the  other  end  was  either  free,  supported,  clamped,  or 
constrained  in  a  non-standard  way  (we  term  this  end  the  “constrained”  end).  An  impulse 
force  was  applied  to  the  stationary  beam  at  the  free  end,  and  the  resulting  deformation 
and  slope  of  the  constrained  end  measured.  Knowing  these  deformations  and  slopes  is 
equivalent  to  knowing  three  sets  of  natural  frequencies,  which  are  necessary  to  infer  the 
elastic  properties  of  the  beam.  Barcilon  refers  to  these  trio  of  spectra  as  “sympathetic” 
spectra.  Given  these  three  sympathetic  spectra,  the  solution  of  the  inverse  problem,  if  it 
exists,  is  unique.  In  a  later  paper,  Barcilon  [4]  investigated  an  apparent  paradox  between 
this  uniqueness  result  and  a  paper  by  Boley  and  Golub  [6],  who  find  a  multiplicity  of 
solutions  when  constructing  a  symmetric  pentadiagonal  matrix  from  its  spectra.  Barcilon 
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states  that  beam  vibration  problems  lead  to  peiitadiagonal  systems  when  recast  as  finite 
difference  problems,  and  hence  the  apparent  paradox.  This  apparent  paradox  is  resolved 
when  Barcilon  shows  that  Boley  and  Golub  chose  their  three  spectra  without  regard  for 
“sympathy”  between  the  spectra. 

Gladwell  [8]  modelled  a  beam  using  rigid  rods  joined  together  by  rotational  springs, 
with  lumped  masses  at  the  joints.  One  end  of  the  beam  was  clamped,  while  the  other 
end  met  one  of  four  conditions:  free,  pinned,  sliding,  or  clamped.  Gladwell  establishes 
necessary  and  sufficient  conditions  for  the  existence  of  a  discrete  model  having  a  given 
spectrum,  and  sets  up  a  procedure  to  find  the  model.  In  a  later  paper  Gladwell  [9]  reviewed 
the  literature  for  solutions  to  inverse  vibration  problems.  Essentially  that  review  looked 
at  the  problem  of  determining  the  system’s  properties  (for  example,  mass  and  stiffness) 
from  vibration  measurements. 

Berman  [5]  investigated  the  problem  of  system  identification  using  data  obtained  from 
dynamic  tests  of  the  structure.  The  structure  was  modelled  by  a  linear  mass,  damping, 
and  stiffness  matrix.  Berman  concludes  that  the  most  promising  approach  to  modelling  is 
to  use  test  data  to  minimally  modify  a  realistic  analytic  model  (subject  to  a  set  of  physical 
constraints). 

6ry,  Glaser,  and  Holzdeppe  [22]  reconstructed  external  and  internal  forces  based  on 
measured  structural  responses.  They  assumed  a  priori  knowledge  of  the  mass  distribution 
and  dynamic  behaviour  of  the  system,  and  a  linear  elastic  system  with  proportional  viscous 
damping.  The  number  of  dynamic  response  measuring  locations  should  be  higher  than 
the  number  of  significant  modes.  They  gave  an  example  of  a  discretised  cantilever  beam 
that  had  several  measurement  locations  along  the  beam.  This  work  was  extended  in  later 
papers  by  one  of  the  above  authors  [23,  24]. 

Arosio,  Panizzi,  and  Paoli  [2]  proved  the  well-posedness  of  a  Timoshenko  beam^  with 
axially  varying  physical  properties  and  sliding  ends.  They  state  that  the  equation  cannot 
be  studied  using  an  iteration  of  the  Fourier  series,  and  instead  used  a  variational  approach 
developed  by  Washizu. 

The  dynamics  of  the  Timoshenko  beam  were  recast  by  Gopalakrishnan,  Martin,  and 
Doyle  [11]  so  that  the  description  only  requires  information  at  the  end  points.  The  re¬ 
sulting  dynamic  stiffness  relations  were  assembled  (akin  to  finite  elements)  allowing  exact 
frequency  dependent  response  for  the  Timoshenko  beam  irrespective  of  element  length. 

Using  a  boundary  integral  equation  Tanaka  and  Bercin  [26]  developed  the  solution 
for  the  free  vibration  of  a  Timoshenko  beam.  A  general  Timonshenko  beam  of  open 
cross-section  with  non-coincident  shear  centre  and  centroid  was  modelled.  They  showed 
that  unacceptably  large  errors  result  from  the  simpler  Bernoulli-Euler  beam  theory  model 
(especially  for  higher  order  modes).  Lee  and  Lin  [18]  developed  an  approximate  solution 
for  the  transverse  vibration  of  a  non-uniform  Bernoulli-Euler  beam  with  time-dependent 
elastic  boundary  conditions. 

The  vibration  characteristics  of  blades  with  multiple-load-paths  (at  the  root)  were 
determined  by  Lauzon  and  Murthy  [17].  They  developed  a  modified  Galerkin’s  method 
to  model  a  non-rotating  beam  undergoing  coupled  fiapwise  bending,  chordwise  bending, 

'  Timoshenko  beam  theory  corrects  for  the  rotary  inertia  and  shear  in  classical  beam  theory  [26] . 
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twisting,  and  extensional  motions.  The  derived  natural  frequencies  of  a  simulated  rotor 
compared  well  with  both  experimental  results  and  a  finite  element  approach. 

A  thorough  survey  and  comparison  of  engineering  beam  theories  for  helicopter  rotor 
blades  was  undertaken  by  Kunz  [16].  The  characteristics  and  differences  amongst  the 
various  formulations  were  reviewed,  and  the  results  presented  in  a  historical  context. 

D’Cruz  [7]  determined  the  location  of  forces  on  a  plate.  First  the  magnitude  and 
location  of  a  static  point  load  on  an  infinite,  and  then  semi-infinite,  rectangular  plate 
were  determined.  This  procedure  was  then  extended  to  harmonic  point  loads.  These  load 
magnitudes  and  locations  were  determined  using  arrays  of  discrete  points  on  the  plate.  The 
solution  involved  the  least-squares  minimisation  of  an  objective  function,  which  was  linear 
in  the  force  magnitude.  This  load  determination  problem  was  found  to  be  ill-conditioned 
near  boundaries. 


2  Equations  Governing  the  Vibration  of  a 

Cantilever  Beam 


In  this  section  we  develop  the  PDF  governing  the  transverse  vibration  of  a  beam. 


Figure  2.1:  A  cantilever  beam  with  a  non-uniform  loading  distribution 
f,  which  varies  with  space  and  time.  The  beam  has  length  I,  mass  per 
unit  length  m,  and  stiffness  El. 


We  begin  by  considering  a  cantilever  beam  of  length  I  with  load  per  unit  length  /(x,  i) 
that  varies  both  spatially  and  temporally  (see  Figure  2.1).  The  spatial  dimensions  x  and  y 
define  the  distance  along  the  beam  and  the  beam’s  transverse  deflection  respectively.  The 
temporal  dimension  is  denoted  by  t.  Non-dimensionalising  these  quantities  we  obtain  the 
non-dimensional  variables 


t  = 


and 


y 

y=v 

P  - 


where  m,  F,  and  I  are  the  beam’s  mass  per  unit  length.  Young’s  modulus,  and  second 
moment  of  area  (or  cross-sectional  area  moment  of  inertia)  respectively.  For  simplicity  we 
assume  the  mass  m  and  stiffness  El  are  constant  along  the  beam. 
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We  see  from  Figure  2.1  that  the  boundary  conditions  at  the  fixed  end  {x  =  0)  are 


y(0,  t)  =  0  and 


dx 


=  0. 


x=0 


(2.1) 


These  two  conditions,  termed  “geometric”  boundary  conditions,  impose  zero  deflection 
and  zero  slope  at  the  fixed  end  of  the  beam.  Similarly  we  see  that  the  free  end  has  the 
boundary  conditions 


d'^y 

dx'^ 


=  0 


and 


.T=l 


d^y 

dx^ 


=  0. 


x=l 


(2.2) 


These  two  conditions,  termed  “natural”  boundary  conditions,  impose  zero  bending  mo¬ 
ment  and  zero  shear  at  the  free  end  of  the  beam.  This  last  statement  follows  from  the  fact 
that  (perpendicular  to  the  beam)  the  bending  moment  and  shear  are  respectively  given  by 
El  d'^y/dx^  and  El  d^y/dx^.  For  further  details  on  the  relation  between  bending  moment 
or  shear  and  derivatives  of  beam  deflection  see,  for  example,  Young  [30]  or  Timoshenko 
and  Gere  [27]. 


We  denote  the  initial  conditions  as 


y{x,0)  =  yo{x)  and 


no(.x), 


(2.3) 


that  is,  yo  and  vq  denote  the  initial  deformation  and  velocity  respectively. 

Using  Meirovitch  [20,  Eq.  (5.43)]  we  know  that  the  PDE  governing  the  transverse 
vibration  of  a  beam  with  constant  mass  and  stiffness  is  given  by 


d'^y  d‘^y 


0  <  X  <  I, 


or  in  non-dimensional  variables 


d'^y  d^y 
dx^  dt'^ 


0  <  a;  <  1. 


(2.4) 


The  following  assumptions  were  made  in  developing  the  PDE  shown  above: 

•  The  rotation  of  the  differential  elements  is  insignificant  (simple  beam  theory). 

•  The  shear  deformation  is  small  in  relation  to  the  bending  deformation  (simple  beam 
theory). 

•  The  ratio  between  the  beam’s  length  and  height  is  relatively  large  (say  more  than  10). 

•  The  beam  does  not  become  too  “wrinkled”  in  flexure. 

Higher  modes  obtained  by  using  Equation  (2.4)  are  inaccurate  because  both  the  simple 
beam  and  the  “wrinkling”  assumptions  are  violated.  One  way  to  determine  higher  mode 
solutions  would  be  to  numerically  solve  the  exact  PDE  governing  beam  vibration  (that  is, 
without  the  above  assumptions). 
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3  Separation  of  Variables  Solution 


In  this  section  we  solve  the  PDE  of  a  vibrating  cantilever  beam,  developed  in  the 
previous  section,  using  the  separation  of  variables  technique.  For  an  explanation  of  this 
technique  see,  for  example,  O’Neil  [21].  We  follow  a  similar  analysis  by  Meirovitch  [20] 
(specifically  pages  223-227  and  235-238),  who  determines  the  modal  response  of  a  hinged 
beam.  (A  hinged  beam  is  one  that  has  pin  joints  at  either  end  and  no  joints  along  the 
beam.) 

We  assume  the  solution  may  be  partitioned  into  a  space  dependent  function  Y  =  Y(x) 
and  a  time  dependent  function  T  =  T(t),  that  is 

y(x,t)  =  Y(x}T(t). 

Substituting  the  above  expression  into  the  homogeneous  form  (that  is,  /  =  0)  of  the 
vibration  PDE  (Equation  (2.4))  gives 

y(4)  J,  +  yf  ^  0, 


and  hence 

y(4)  y 


where  /?  is  a  constant  from  the  separation  of  variables  technique.  The  over-dot  denotes 
differentiation  with  respect  time  (for  example,  T  =  (fT/dt^);  while  the  dash  and  bracketed 
superscript  denote  differentiation  with  respect  to  space  (for  example,  Y”  =  d?Y / dx"^  and 
=  d^YJdx^).  The  constant  of  separation  /?  is  related  to  the  natural  frequency  of  the 
dimensional  beam  y/ E I / (see,  for  example,  Meirovitch  [20,  p.  226]). 

The  solution  of  the  fourth  order  differential  equation  Y^^"^  —  (3'^  Y  =  0  is  given  by 
Tuma  [28,  p.  181]  as 

Y{x)  —  A cosh(/9a:)  -|-  B  sinh(/?x)  -f  C cos{/3x)  +  D  sm{f3x), 


where  A,  B,  C,  and  D  are  constants  to  be  determined  from  the  boundary  conditions. 
Using  the  boundary  conditions  given  by  Equations  (2.1)  and  (2.2)  we  have  that 


y(0)  =  0,  y'(o)  =  0, 

y"(l)=0,  and  y(3)(l)  =  0. 

Solving  for  the  four  constants  using  the  above  boundary  conditions  yields 
Yj{x)  =  Aj  [(cosh  PjX  -  cos  Pjx)  —  aj  (sinh  PjX  —  sin  Pjx)] , 
where  is  the  jth  positive  root  of  the  equation 


cos  Pj  cosh  I3j  =  —  1 , 


and 


cosh  Pj  +  cos  pj 
sinh  Pj  +  sin  Pj 


(3.1) 


(3.2) 


(3.3) 


is  a  constant,  which  we  term  the  hyperbolic  constant.  In  the  above  solution  we  have 
ignored  all  answers  leading  to  the  trivial  solution  Y{x)  =  0  (such  as  Aj  =  0).  Note  that 
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all  the  solutions  Yj  are  orthogonal,  that  is,  Yi{x)Yj{x)dx  =  0  ifi  ^  j.  Kelly  [14,  p.  478] 
and  Meirovitch  [19,  p.  163]  give  the  same  expression  for  the  solution  of  a  cantilever  beam 
under  free  vibration  as  Equation  (3.2). 

We  want  to  select  Aj  so  that  the  solutions  are  not  only  orthogonal  but  are  in  fact 
orthonormal,  that  is, 

[  Yi{x)Yj{x)dx^\^  (3-4) 

Jo  [1  ifi=;. 

Using  Mathematica  [29]  (a  computer  program  that  can  perform  symbolic  manipulation) 
the  constants  Aj  that  satisfy  this  orthonormal  condition  were  determined.  Some  tedious 
algebraic  manipulation  (using  Mathematica  and  Abramowitz  and  Stegun  [1])  then  led  to 
the  solution 

^2 _ iJj  (sinh/?j  +  sixiPj)"^ _ 

^  ~  (dj  (sinh (3j  +  sin  pjf  +  3  (1  +  cos (3j  cosh  f3j)  (cos Pj  sinh  pj  -  cosh  pj  sin  Pj) 

Using  the  first  part  of  Equation  (3.3)  the  constants  Aj  simplify  to 

Aj  =  l.  (3.5) 


Thus  Ec^uation  (3.2)  becomes 

Yj{x)  =  (cosh PjX  -  cos pjx)  -  aj  (sinh  pjX  -  sin Pjx) .  (3.6) 

The  first  four  modes  (using  the  above  equation)  are  shown  in  Figure  3.1. 


X 


Figure  3.1:  The  first  four  vibration  (or  natural)  modes  of  a  cantilever 
beam  (that  is,  Yj{x)  for  j  =  1,2, 3, 4j. 


Before  concluding  this  section  let’s  investigate  the  solutions  of  the  equation  governing 
the  separation  constant  (the  first  part  of  Equation  (3.3)).  Re-arranging  this  equation  gives 


cos  Pj  — 


cosh  Pj 
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and  we  now  clearly  see  that  for  f3j  '>  1  this  equation  is  well  approximated  by  cos/3j  =  0. 

Expanding  the  above  displayed  equation  in  a  Taylor’s  series  about  /3j  =  ttU  — 
ignoring  all  terms  of  O  [Pj  -  'K{j  —  and  above,  and  solving  for  Pj  yields  the  excellent 
approximation 

-  I)  +  2(-l)^sech  [7r(j  -  ^)]  (3.7) 

^  7r(j  -  5)  -  2(-l)^  exp  [-7r(j  -  i)]  (3.8) 

~  Aj  -  5))  (3-9) 

where  the  hat  notation  denotes  an  approximation,  that  is,  Pj  ~  pj. 

The  second  approximation  for  the  constant  Pj  shown  above  comes  from  considering 
the  hyperbolic  function  “sech”  in  exponential  form.  The  symbols  O  defines  the  order,  so 
for  example  if  /(x)  =  0{x‘^)  then  |/(x)/x^|  is  bounded. 

Table  3.1  shows  the  accuracy  of  these  three  approximations  for  Pj  to  six  significant 
figures.  For  each  approximation  both  the  value  and  relative  error  of  the  solution  are  shown, 
where  the  relative  error  is  defined  as  {Pj  —  Pj)/ Pj- 


j 

1 

1.87510 

1.96933 

5.03x10-2 

1.98656 

5.94x10-2 

1.57080 

1.62x10-1 

2 

4.69409 

4.69442 

7.09x10-^ 

4.69440 

7.06x10-^ 

4.71239 

3.90x10-® 

3 

7.85476 

7.85476 

7.66x10-® 

7.85476 

7.67x10-® 

7.85398 

9.88x10“® 

4 

10.9955 

10.9955 

1.02x10-^° 

10.9955 

1.02x10-^° 

10.9956 

3.05x10-® 

5 

14.1372 

14.1372 

1.49x10-^^ 

14.1372 

1.49x10-1® 

14.1372 

1.03x10-^ 

Table  3.1:  Values  and  errors  of  three  approximations  for  the  separation 
constant  Pj .  The  superscripts  ”,  '“t  ”,  and  ”  represent  the  approxi¬ 
mations  given  by  Equations  (3.7),  (3.8),  and  (3.9)  respectively. 


Using  a  similar  approach  for  the  constant  cxj  yields  the  excellent  approximation 

5,^ _ 1 _  2(-l)4ech  [7T{j  -  I)] 

tanh  [7r{j  -  i)]  -  (-l):'sech  [Tr{j  -  i)]  +  sinh  [7r(j  -  ^)] 

1  +  2(-l)-^exp  [-7r(ji  -  ^)]  (3.11) 

«1.  (3.12) 


In  Equation  (3.10)  we  have  made  use  of  the  approximation  for  the  separation  constant  pj 
given  by  Equation  (3.7).  Table  3.2  shows  the  values  and  relative  errors,  (Sj  —  aj)/aj,  of 
these  three  approximations  to  six  significant  figures. 
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3 

s; 

\aj  -  S;|/Oj 

j 

\aj  -  S/  /oj 

0 

IM 

1 

0.73410 

0.68692 

6.43x10-2 

0.58424 

2.04x10-1 

1.00000 

3.62x10 

-1 

2 

1.01850 

1.01850 

5.76x10-® 

1.01800 

4.92x10-^ 

1.00000 

1.81x10 

-2 

3 

0.99922 

0.99922 

4.68x10-1® 

0.99922 

9.03x10-^ 

1.00000 

7.76x10 

-4 

4 

1.00000 

1.00000 

3.77x10-1^ 

1.00000 

1.69x10-® 

1.00000 

3.36x10 

-5 

5 

1.00000 

1.00000 

3.05x10-1® 

1.00000 

3.15x10-12 

1.00000 

1.45x10 

-6 

Table  3.2:  Values  and  errors  of  three  approximations  for  the  hyperbolic 
constant  aj.  The  superscripts  ""j-”,  and  represent  the  approxi¬ 
mations  given  by  Equations  (3.10),  (3.11),  and  (3.12)  respectively. 


As  can  be  seen  from  Tables  3.1  and  3.2  the  approximations  for  the  separation  and 
hyperbolic  constants  are  excellent  (especially  for  j  >  2). 


4  Solution  in  terms  of  Different 
Loading  Distributions 

We  now  consider  several  different  loading  conditions.  We  will  begin  by  looking  at  the 
simplest  loading  (impulse  loading  in  §  4.1),  then  investigate  the  solution  under  a  slightly 
more  complex  loading  pattern  (step  loading  in  §  4.2).  Finally,  in  §  4.3,  we  develop  the 
solution  in  a  more  intrinsic  way  by  making  use  of  the  beam’s  natural  modes. 

Let  us  consider  the  non- homogeneous  solution  to  the  vibrating  beam  problem.  Using  an 
approach  analogous  to  Fourier  series  decomposition  (or  any  orthogonal  decomposition), 
we  may  represent  the  vibration  solution  as  a  series  expansion  of  orthogonal  functions. 
Using  the  homogeneous  (and  orthogonal)  solutions  given  by  Equation  (3.6)  we  can  write 
the  vibration  solution  as 

C» 

i=i 

where  ijjj  are  the  coefficients  in  the  series  expansion  of  the  solution  in  terms  of  vibra¬ 
tional  modes.  We  term  the  coefficients  i/jj  the  solution  coefficients.  Similarly  the  loading 
distribution  has  the  series  expansion 

00 

f{x,t)  =  Y^(t)j{t)Yj{x),  (4.2) 

j=i 

where  (t>j{t)  is  the  coefficient  for  the  jth  homogeneous  solution.  We  term  the  coefficients 
(j)j  the  loading  coefficients.  Note  that  non-zero  derivatives  at  the  fixed  end,  that  is 

/  0,  for  *  =  0, 1, 2, . . ., 

x=Q 


dx^ 
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are  not  efficiently  represented  by  the  formulation  of  Equation  (4.2)  (see  Figure  3.1).  In 
other  words,  even  simple  functions  such  as  cos  x  would  need  a  large  number  of  modes  to 
be  well  approximated. 

If  we  substitute  the  vibration  solution  given  by  Equation  (4.1)  and  the  loading  dis¬ 
tribution  given  by  Equation  (4.2)  into  the  PDE  governing  beam  vibration  (2.4)  we  have 
that 

00  OO  OO 

j  =  l  j=l  j  =  l 

Since  the  fourth  derivative  of  the  jth  mode  is  a  multiple  of  the  mode  itself  (that  is, 

=  (3^  Yj),  we  have  that 

OO 

+  /?  -  ^j)  =  0- 

3  =  3 

The  modes  Yj  are  orthogonal,  however,  so  we  must  have  that 

i’j  (t)  +  (if  ipj  (t)  =  ((j  (t) ,  (4.3) 

for  all  j.  Solving  this  second  order  differential  equation  (DE)  for  the  solution  coefficients 
'tpj  determines  the  vibration  response  of  the  cantilever  beam. 

4.1  Non-Uniform  Impulse  Loading  Along  the  Beam 

Consider  the  discrete  loading  distribution  shown  in  Figure  4.1.  This  loading  configu¬ 
ration  may  be  represented  mathematically  as  a  set  of  spatial-impulse  loads,  that  is, 

K 

=  (4.4) 

i=i\ 

where  and  are  the  magnitude  and  location  of  the  ith  discrete  load  respectively,  K  is 
the  number  of  discrete  spatial- impulse  loads,  and  5{a)  is  the  delta  (or  impulse)  function 
defined  by  S{x  —  a)f{x)dx  =  f{a)  for  any  e  >  0.  Note  that  both  the  loadings  and 
locations  have  been  appropriately  non-dimensionalised  (see  page  3),  and  that  the  load 
magnitudes  vary  temporally.  Figure  4.1  shows  these  new  definitions  superimposed  on 
our  cantilever  beam.  Finally,  note  that  Equation  (4.4)  is  composed  of  a  set  of  impulse 
functions  5{x  —  ^i)  only  in  space  and  not  in  time. 

Multiplying  both  sides  of  Equation  (4.2)  by  Yk{x)  and  integrating  with  respect  to  x 
over  the  length  of  the  beam  gives  that 

/  y]<Pj(.tWj{^)yk{x)dx=  /  f{x,t)Yk{x)dx, 

Jo  ^  Jo 

where  we  have  switched  the  left  and  right  hand  sides  as  compared  to  Equation  (4.2). 
Assuming  the  series  given  by  Equation  (4.2)  is  “well  behaved”  we  can  change  the  order  of 
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Figure  4-1-  cantilever  beam  with  spatial-impulse  loading  $i. 


summation  and  integration.  (For  a  more  detailed  discussion  on  the  integration  of  series 
see  Knopp  [15,  p.  341].)  Furthermore,  using  Equation  (4.4)  in  the  right  hand  side  yields 

oo  /-l  pi  ^ 

V (Pj{t)  /  Yj{x)Yk{x)dx  =  /  y^S{xi-^i)^i{t)Yk{x)dx 

Jo  Jo 

oo  K  pi 

Sjk  =  Y  j  yk{x)dx, 

j=i  j=i  ° 

where  Sjk  is  the  Kronecker  delta  function,  and  by  definition  Skk  =  1  and  Sjk  =  0  if  j  7^  fc. 
In  the  left  hand  side  of  the  above  equations  we  have  made  use  of  the  orthonormal  property 
of  the  homogeneous  solution  (see  Equation  (3.4))  to  get  from  the  first  line  to  the  second 
line.  We  can  change  the  order  of  integration  and  summation,  on  the  right  hand  side,  since 
both  operations  are  finite.  Using  the  definition  of  the  delta  generalised  function  we  have 
that 

K 

0,(t)- j]$,(t)n.(6).  (4.5) 

i=l 

Thus  our  spatial-impulse  loading  (4.4)  can  be  written  as  the  sum  of  the  beam’s  vibrational 
modes. 

The  analysis  in  the  remainder  of  this  subsection  can  be  greatly  simplified  if  we  assume 
the  loading  is  constant  in  time;  we  will  later  explain  the  need  for  this  restriction.  By 
“constant  in  time”  we  mean  that  the  loading  is  applied  as  a  step  load,  so  that  the  loading 
is  constant  both  before  and  after  its  application. 

Under  a  constant  temporal  loading  assumption,  substituting  Equation  (4.5)  into  Equa¬ 
tion  (4.3),  and  solving  for  the  resulting  second  order  DE  gives  the  solution 

1  ^ 
o'  i=l 

where  Aj  and  Bj  are  the  constants  from  the  solution  of  the  DE. 
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Using  the  initial  conditions  '^j(O)  =  0  and  '0j(O)  =  0  (that  is,  zero  initial  deformation 
and  velocity)  we  have  that 

1  ^ 

n  i=i 

Hence  the  vertical  deformation  of  a  vibrating  beam  with  zero  initial  conditions  and  spatial- 
impulse  loading  is  given  by 

oo  ^  K 

vM  =  E  ^  [1  -  ““W"  ‘)1  E  (6) 

j=l  n  i=l 

We  can  now  begin  to  see  that  having  vibration  information  at  one  point  on  the  beam 
(say  from  an  accelerometer)  would  allow  us  to  determine  the  loading  distribution  along 
the  beam.  To  see  this  more  clearly  differentiate  the  solution  given  by  Equation  (4.6)  twice 
with  respect  to  time,  which  gives  the  acceleration  of  the  beam  as 

K  oo 

i=l  j=l 

Changing  the  order  of  summation  in  the  above  expression  may  be  technically  incorrect 
(using  classical  definitions  of  convergence),  since  it  is  not  at  all  clear  that  the  series  is 
convergent.  In  fact  (using  classical  definitions  of  convergence)  the  series  may  prove  to  be 
divergent,  since  for  j  large  we  have  that  Yj{x)  ^  sin[7r(_)  —  ^)x]  —  cos[7r(j  —  ^)2;].  (For 
an  interesting  discussion  on  divergent  series  see  Hardy  [13],  who  states  that  depending  on 
definitions  some  “classically”  divergent  series  are  in  fact  convergent.)  As  an  aside  note 
that  this  “non-convergence”  result  is  not  unexpected,  remember  that  an  impulse  function 
excites  all  frequencies  of  a  Fourier  series  equally.  The  fact  that  the  impulse  function  is  a 
generalised  function  (see  for  example  Greenberg  [12])  should  alert  us  to  these  apparently 
unusual  results. 

If  (for  the  sake  of  the  argument)  we  assume  the  series  to  be  convergent  on  physical 
considerations,  then  the  above  re-arrangement  is  permissible  and  the  equation  displayed 
above  can  be  re-written  in  the  form 

K 

y{x,t)  =  J^7i(a;,t)$i,  (4.7) 

i=l 

where 

OO 

t)  =  ^  cos(/3/ 1)  Yj{x)  Yj{^i). 
i=i 

An  example  will  crystalise  the  idea  of  determining  the  magnitude  of  the  loads  3>i.  Let’s 
take  K  samples  of  acceleration,  at  times  t  ~  ,  tx,  from  an  accelerometer  located 

at  X  =  along  the  beam.  As  an  example  let  K  —  3,  then  from  Equation  (4.7)  we  have 
the  linear  system  of  K  equations 


’71(6,  ii)  72(C*,^i)  73(6,  ti)’ 

r$i 

y{^*,i2) 

= 

71  ('^*5^2)  72(^*,i2)  73(^*,^2) 

$2 

7l(^*,^3)  'y2{^*,h) 

$3 
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Since  the  vector  of  accelerations  on  the  left  and  the  matrix  of  constants  on  the  right  are 
both  known,  the  above  linear  system  can  easily  be  solved  for  the  vector  of  unknown  load 
magnitudes  on  the  right. 

If  the  location  of  the  applied  loads  were  also  unknown  then  Equation  (4.7)  would  yield 
a  non-linear  system  of  equations  in  terms  of  and  ^i.  As  can  be  seen  this  non-linear 
system  of  equations  would  require  2K  acceleration  samples  to  solve  for  the  magnitude  and 
location  of  the  unknown  loads.  Moreover,  because  the  system  is  non-linear  there  is  no 
guarantee  that  the  solution  would  be  unique. 

Consider  the  following  thought  experiment.  First  uniformly  space  the  loading  locations 
along  the  beam,  and  let  the  number  of  these  locations  K  tend  to  infinity.  Secondly,  using 
an  accelerometer  record  the  vibration  at  one  point  along  the  beam.  We  can  then  determine 
the  loading  distribution  along  the  entire  beam  using  vibration  information  from  only  one 
point  on  the  beam!  This  thought  experiment  assumes  we  can  record  all  frequencies  (that 
is,  an  infinite  number  of  natural  freciuencies) .  In  reality  we  would  be  restricted  by  the 
bandwidth  capability  of  our  measuring  system,  which  would  in  turn  restrict  the  accuracy 
of  our  load  reconstruction  capability.  Remember  also  that  we  have  made  the  assumption 
that  the  beam  doesn’t  become  too  “wrinkled” ,  and  hence  we  could  not  accurately  model 
higher  modes. 


4.2  Non-Uniform  Step  Loading  Along  the  Beam 


Consider  the  discrete  loading  distribution  shown  in  Figure  4.2.  This  loading  configu¬ 
ration  may  be  represented  mathematically  as  a  set  of  hat  loads,  that  is. 


K 

f{x,  =  6, 6+i)'*’i(0i 

i=0 


where  for  alH  =  0, 1, . . . ,  ii',  and 


/i(x;  a,  b) 


0  if  a;  <  a, 

<  1  if  a  <  X  <  6, 
0  if  X  >  6, 


(4.8) 


(4.9) 


is  the  hat  function.  The  hat  function  loading  representation  given  by  Equation  (4.8)  uses 
the  two  spatial  locations  and  which  we  define  as  x  =  0  and  x  =  1  (the  fixed  and 

free  ends)  respectively.  Again  both  the  load  magnitudes  and  locations  have  been 
appropriately  non-dimensionalised . 

Following  a  procedure  that  is  analogous  to  the  one  used  in  the  previous  subsection  we 
arrive  at  the  solution  for  a  hat  loaded  beam.  For  zero  initial  deformation  and  velocity 
(yg  =  uo  =  0)  the  vertical  deformation  of  a  cantilever  beam  is  given  by 

oo  .  K 

yM  =  [1  -  (4-10) 

j=\  n  i=0 
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Figure  J^.2:  A  cantilever  beam  with  step  loading 


where 


[  (sinh  PjX  —  sin  Pjx)  —  aj  (cosh  f3jX  +  cos  f^x)  ] 


?i+l 


is  the  integral  of  the  jth  natural  mode  over  the  ith  beam  segment. 


4.3  Solution  in  Terms  of  Natural  Modes 

The  problem  of  choosing  beam  segmentations  may  be  alleviated  by  deriving  beam 
loading  results  in  terms  of  the  beam’s  natural  frequencies.  In  this  subsection  we  reconstruct 
the  loading  along  a  beam  in  terms  of  the  beam’s  modes.  (We  will  see  that  this  type  of 
construction  leads  to  a  more  natural  solution.)  Figure  4.3  shows  an  example  of  this  natural 
mode  loading  for  the  distribution  1.67Fi(a;)  +  0.190y2(x)  +  0.0889y3(x)  +  0.390l4(x).  The 
complete  load  is  shown  as  the  dark  grey  region,  while  the  individual  decomposed  modes 
are  shown  as  dashed  lines. 

In  an  analogous  way  to  the  previous  two  subsections  we  can  express  the  loading  coeffi¬ 
cients  (bj  in  terms  of  the  loading  function  /(x,  t).  Multiplying  both  sides  of  Equation  (4.2) 
by  a  natural  mode  and  integrating  over  the  length  of  the  beam  gives  the  jth  loading 
coefficient  as  ^ 

=  f  f{x,t)Yj{x)dx.  (4.11) 

Jo 

Assume  the  loading  is  periodic  in  time  with  period  P,  that  is,  f{x,t  +  P)  =  f{x,t). 
The  loading  coefficient  may  then  be  expanded  in  a  Fourier  series 

00 

(t)j{t)  =  aoj  +  y]  [atj  cos{uJit)  +  bij  sin(c<;it)] ,  (4-12) 
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Ooi 


^^03 

002 

O04 


Figure  4-3:  A  cantilever  beam,  with  the  loading  decomposed  into  the 
beam.’s  natural  modes.  The  dashed  lines  represent  the  individual  de¬ 
composed  modes,  while  the  dark  grey  region  represent  the  sum  of  these 
m,odes  (that  is,  the  total  loading). 


where  Oy  and  bij  are  the  Fourier  coefficients  for  the  cosine  and  sine  expansions  respectively 
and  oJi  =  {2Tri)/P  is  the  ith  loading  frequency. 

The  solution  to  the  second  order  differential  equation  given  by  Equation  (4.3)  with  the 
above  equation  substituted  for  the  loading  coefficient  is  given  by 


^  Aj  cos(/3/ 1)  -f  Bj  sin(/?/  0  +  ^ 


aij  cos(a;jt)  -f-  bij  sin(wit) 


(4.13) 


(see,  for  example,  Tuma  [28,  p.  180]). 

Once  again  using  the  orthonormal  property  of  the  normal  modes  we  have  from  Equa¬ 
tion  (4.1)  that 

/  y{x,t)Yj{x)dx. 
do 

Hence  for  the  initial  condition  (t  —  0)  we  have  equating  the  above  expression  with  Equa¬ 
tion  (4.13) 

OO  ^ 

"  "if  " 

where  ^ 

aj  =  /  y{x,0)Yj{x)dx  and 
do 

are  respectively  the  coefficients  of  the  beam’s  initial  deformation  and  velocity  decomposed 
into  the  beam’s  natural  modes.  Substituting  the  above  expressions  for  the  constants  Aj 
and  Bj  into  Equation  (4.13)  gives  the  final  form  of  the  series  expansion  coefficients 


r,- 


-E 


Uj 


=  [  yix,0)Yj{x)dx 
do 


Vn(t) 


cos(/i^  t) 


CT  J 


f) 


+  E 

1=1 


cos(tujt)  —  cos{l3^t) 

Aj  -  A) 


Ti  + 


an  + 


—  {1  -  cosiBf  t)}  \  aoj 
sin(u;jt)  -  {ui / (df)  sm{j5h) 

W^) 


.  (4.14) 
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We  want  to  determine  how  the  accelerations  at  one  location  along  the  beam  relate  to 
the  beam  loading.  Differentiate  Equation  (4.1)  twice  with  respect  to  time  to  yields 

OO 

1=1 

which  shows  that  the  acceleration  of  the  beam  is  related  to  the  acceleration  of  the  solu¬ 
tion  coefficients  V-’i-  Differentiating  Equation  (4.14)  twice  with  respect  to  time  gives  the 
acceleration  of  the  solution  coefficients  as 


cos{l3^t)]  aj  +  [-13^  tj  +  [cos{(3^t)]  aoj 


+E 

i=l 


cos{Pjt)  —  ojf  cos{u>it) 


(/? 


UJ. 


ujil3f  sm{l3h)  -  Jf  sin(cuit) 


00. 


'■) 


(4.16) 


In  the  above  expression  all  the  quantities  in  square  brackets  are  known.  At  first  it  would 
appear  we  can  reconstruct  the  loading  along  the  beam  if  we’re  given  a  set  of  acceleration 
measurements  at  one  location  on  the  beam.  Unfortunately,  the  constant  loading  (that 
is,  the  coefficients  aoj)  and  the  initial  deformation  (crj)  are  entangled,  and  using  a  single 
measurement  location  these  unknowns  cannot  be  disentangled. 

The  simplest  way  to  see  this  entanglement  between  the  initial  deformation  Oj  and  the 
constant  force  aoj  is  to  apply  two  substitutions.  Making  the  substitutions 

aj  -f  Act  and  aoj  +  f3j  Aa 

for  aj  and  aoj  respectively  in  Equation  (4.15)  and  simplifying  returns  exactly  the  same 
expression  (namely,  Equation  (4.15))  regardless  of  the  value  of  Act.  The  fact  that  two 
different  initial-deformation  and  temporally-constant-loading  combinations  give  the  same 
acceleration  for  the  solution  coefficient,  means  that  these  two  unknowns  {aj  and  aoj) 
cannot  be  disentangled  without  additional  information. 

Would  using  strain  instead  of  acceleration  improve  this  entangled  situation?  Yes — but 
only  mildly,  since  we  have  only  gained  one  piece  of  additional  information  instead  of  the 
n  pieces  of  information  we  require  to  approximate  the  loading  using  the  first  n  modes. 

We  know  that  the  strain  is  proportional  to  the  bending  moment,  and  hence  strains  are 
also  proportional  to  d'^y/dx^.  (The  constant  of  proportionality  between  strain  and  this 
second  derivative  is  the  distance,  measured  in  the  bending  plane,  from  the  beam’s  neutral 
axis  to  the  desired  strain  location.)  Differentiating  Equation  (4.1)  twice  with  respect  to 
the  horizontal  distance  along  the  beam  we  have  that 


dx^ 


'^i)j{t)Yj"{x), 

t=i 


(4.16) 


and  from  Equation  (3.6)  the  second  derivative  of  the  jth  natural  mode  is 


Yj"{x)  =  [cosh  {f3jx)  —  cos  {Pjx)]  —  ocj  [sinh  {(3jx)  —  sin  {Pjx)]  |. 
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If  we  use  strain  instead  of  acceleration  to  determine  the  loading,  then  we  would  use 
the  solution  coefficients  themselves  (Equation  (4.14))  instead  of  their  accelerations  (Equa¬ 
tion  (4.15)).  As  already  mentioned  above,  using  strain  instead  of  acceleration  improves  the 
situation  slightly  (it  gives  us  one  additional  piece  of  information).  To  see  that  the  initial 
deformation  and  temporally-constant  force  are  still  entangled  we  will  again  use  substitu¬ 
tions,  but  this  time  use  multiple  modes  instead  of  just  one.  Making  the  substitutions 

(ji  -t-  Act,  (Zoi  +  ~  Y"{^  )^ 

for  (Ti,  floi,  (T'i,  and  aoj  (for  j  >  2)  respectively  in  Equation  (4.16)  and  simplifying  returns 
exactly  the  same  expression  (namely.  Equation  (4.16))  regardless  of  the  value  of  Act.  As 
defined  earlier,  the  spatial  position  represents  the  location  of  the  measurement  device 
(in  this  case  a  strain  gauge). 

Near  the  fixed  end  of  the  beam,  from  an  implementation  point  of  view,  strain  mea¬ 
surements  would  be  preferable  to  acceleration  measurements,  since  strains  would  tend  to 
be  high  in  this  region  while  accelerations  would  be  low.^ 

In  a  similar  way  we  can  show  that  the  Fourier  coefficients  Oy  and  fey  are  entangled 
with  the  initial  deformation  coefficient  aj  and  initial  velocity  coefficient  Tj  respectively. 

In  summary,  we  have  shown  that  using  acceleration  or  strain  at  one  location  alone  is 
not  enough  to  simultaneously  determine  the  initial  conditions  and  loading  distribution. 
In  fact,  we  cannot  determine  loading  distributions  that  vary  with  time — even  if  we  know 
the  initial  conditions.  We  can,  however,  determine  the  temporally-constant  loading  dis¬ 
tribution  along  the  beam  provided  we  know  the  initial  deformation  (note  that  we  do  not 
require  the  initial  velocity  to  determine  this  loading  distribution). 


5  Simulation 


In  this  section  we  report  on  a  simulation  of  load  prediction  on  a  simple  cantilever  beam. 
The  transient  response  of  a  finite  element  (FE)  cantilever  model  is  used  to  determine  the 
beam’s  strain  at  one  location.  To  test  the  load  prediction  technique  developed  in  §  4.3, 
we  predict  the  load  distribution  given  the  transient  strain  response  at  one  location  along 
the  FE  model. 

The  beam  arbitrarily  chosen  to  verify  this  load  prediction  capability  was  a  tube 
1000  ram  in  length,  with  a  wall  thickness  of  10  mm  and  a  square  cross-section  of  50  mm 
by  50  mm  (outer  dimensions).  The  Young’s  modulus  and  density  of  this  tube  were 
72.40  X  10®  N/m^  and  2.768  x  10^  kg/m^  respectively.  (The  Poisson’s  ratio  and  the 
shear  modulus  were  0.3  and  27.94  x  10®  N/m^  respectively.  However,  these  two  properties 
were  not  necessary  for  the  vibration  response  solution.)  The  constant  load  —  x)  was 
applied  as  a  step  function,  where  x  is  the  distance  (in  metres)  along  the  tube  from  the 
fixed  end.  This  load  represents  a  cubic  distribution  with  zero  load  at  both  ends  and  zero 
slope  at  the  fixed  end  of  the  tube.  This  step  load  was  instantaneously  applied  at  fe  =  0, 

^Thanks  to  Shane  Dunn  (a  colleague  at  AMRL)  for  making  me  aware  of  this  point. 
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that  is,  the  load  may  be  represented  as 

x^{l  —  x)  h(t,0,oo),  (5.1) 

where  h{t,0,oo)  is  the  hat  function  defined  in  Equation  (4.9). 

For  simplicity  we  set  the  initial  conditions  to  zero,  that  is,  the  beam  initially  had  zero 
deformation  and  zero  velocity. 

We  will  first  describe  how  we  derived  the  so  called  “exact”  solution  for  this  verification 
problem  in  §  5.1.  We  then  describe  how  the  FE  solution  was  obtained  in  §  5.2.  Finally,  we 
describe  how  successful  the  load  prediction  technique  was,  and  also  discuss  some  practical 
problems  in  §  5.3. 

5.1  The  Exact  Solution 

The  exact  solution  of  a  cantilever  beam’s  response  to  the  step  load  given  by  Equa¬ 
tion  (5.1)  was  found  using  the  solution  procedures  developed  in  §  4.3.  The  magnitudes  of 
the  first  twenty  load  coefficients  {aoj  in  Equation  (4.12))  of  the  exact  solution  are  shown 
in  Figure  5.1. 


mode  number 

Figure  5.1:  The  magnitude  of  the  exact  solution’s  coefficients  (denoted 
by  •)  plotted  against  the  mode  number.  The  relative  difference  ( denoted 
by  0)  between  using  the  corresponding  number  of  modes  as  compared  to 
the  20  mode  solution. 


These  load  coefficients  were  calculated  using  one-hundred  digit  arithmetics.  This  high 
accuracy  was  required  because  of  the  precision  loss  for  the  higher  order  modes  was  signif¬ 
icant.  For  example,  in  calculating  the  load  coefficient  for  the  twentieth  mode  at  least  26 
digits  of  accuracy  were  lost.  The  loss  of  accuracy  is  due  to  the  subtraction  of  numbers 
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with  almost  identical  values,  which  we  see  when  we  re-write  Ecjuation  (3.6)  (the  expression 
for  the  beam  modes)  as 


Yj{x)  =  (1  —  Ofj)  sinh /3jX  -|-  (sin/5jX  —  cosPjx)  -h 


e  —  (1  —  ttj)  sin  PjX 


Despite  the  fact  that  (1  -  aj)  see  Equation  (3.11),  the  first  term 

in  the  above  equation  is  not  necessarily  small  for  j  large.  For  example,  setting  x  =  1  and 
using  the  approximations  given  by  Equations  (3.9)  and  (3.11)  gives  the  approximation 
(1  —  Oj)  sinh  l3jX  «  (— l)-’'*'^  [1  —  which  is  of  order  unity  for  j  large.  This  sort  of 

cancellation  leads  to  the  large  accuracy  losses  in  calculations  involving  higher  modes. 


Also  shown  in  the  log-linear  plot  of  Figure  5.1  is  the  relative  difference  between  using 
only  the  first  n  modes  as  compared  to  using  the  first  twenty  modes.  This  relative  difference 
was  calculated  as 


relative  difference  = 


^fin 


ifin 


ti=0 


\ 


T — [2^20 (^*,  ti  At)] ■ 
^  u=o 


(5.2) 


where  tfin  is  the  number  of  temporal  samples  taken  for  the  evaluation  of  the  relative 
error  and  ynix,t)  are  the  solutions  with  20  modes  and  n  modes  respectively, 

and  At  is  the  time  step  between  temporal  samples.  (In  our  case  t^n  ==  1600  and  the 
non-dimensional  time  step  was  At  =  y/ £T/(m/^)/1600,  see  §  5.2  for  an  explanation  of 
these  value.)  The  relative  difference  defined  by  Equation  (5.2)  is  merely  the  standard 
deviation  of  the  difference  between  the  twenty  mode  and  the  n  mode  solutions,  relative  to 
the  standard  deviation  of  the  twenty  mode  solution.  (This  measure  can  also  be  thought 
of  as  a  root  mean  squared  or  2-norm  measure.)  In  Figure  5.1  the  relative  difference  of 
using  fewer  than  twenty  modes  was  plotted  only  up  to  ten  modes,  which  was  the  number 
of  modes  used  in  the  exact  solution.  In  other  words,  in  the  following  sections  what  we 
term  the  “exact”  solution  contains  only  the  first  ten  modes. 


5.2  The  Finite  Element  (FE)  Solution 

The  FE  package  “Nastran”  was  used  to  solve  for  the  vibration  response  of  our  square 
tubular  beam.^  The  final  FE  model  used  1600  beam  elements  (of  type  “Bar2”),  and  the 
beam’s  response  was  calculated  with  time  steps  of  l/1600th  of  a  second.  The  FE  model 
was  fixed  both  from  moving  in  translation  and  rotation  at  one  end  of  the  beam,  and  was 
free  to  move  at  the  other  end. 

The  loading  was  applied  as  point  loads  at  the  nodes  joining  the  beam  elements,  and 
hence  the  actual  loading  distribution  along  the  beam  can  be  represented  by  a  sum  of 
spatial-impulse  loads.  In  other  words,  the  load  may  be  written  as 

K 

-  x)S{x  -  (i)/K, 

i=0 

^Thanks  to  Soon- Aik  Gan  (a  colleague  at  AMRL)  for  developing  the  Nastran  model. 
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where  is  the  distance  (in  metres)  from  the  fixed  end  to  the  fth  node  and  if  +  1  is 
the  number  of  nodes  {K  =  1600  in  our  case).  The  factor  of  1/if  is  required  in  this  load 
distribution  so  that  the  overall  force  on  the  beam  is  the  same  as  its  continuous  counterpart 
Jq  —  x)dx.  By  integrating  the  above  node  loading  over  the  beam  we  see  that  the  error 
in  the  discretised  node  loading  is  l/(12if^),  and  hence  for  our  1600  element  FE  model  the 
total  load  on  the  beam  has  a  relative  error  of  3.5  x  10~^. 

The  final  time  step  and  number  of  elements  were  decided  upon  after  comparing  sever;, 
solutions  with  varying  time  steps  and  beam  elements.  The  relative  differences  between 
the  solution  with  1600  beam  elements  (and  1600  time  steps)  and  several  other  solutions 
are  shown  in  Table  5.1.  These  relative  differences  were  calculated  as  the  standard  devia¬ 
tion  of  the  difference  divided  by  the  standard  deviation  of  the  1600  element  solution,  see 
Equation  (5.2). 


NUMBER  OF  ELEMENTS 

100 

200 

400 

800 

TIME  STEPS  (sec) 

1/100 

1/200 

1/400 

1/800 

RELATIVE  DIFFERENCE 

9.4  X  10-^ 

9.3  X  10-^ 

2.3  X  10-4 

2.9  X  10~5 

Table  5.1:  The  relative  difference  between  several  different  FE  solutions, 
as  compared  to  the  FE  solution  with  1600  elements  and  a  time  step  of 
1/ 1600th  of  a  second. 


In  order  to  simulate  data  from  a  strain  gauge,  the  stress  on  the  upper  side  of  the  beam, 
50  mm  along  the  horizontal  from  the  fixed  end,  was  extracted  from  the  FE  results.  The 
stress  was  then  scaled  to  the  non-dimensional  strain  result  given  by  the  exact  solution. 
Figure  5.2  shows,  at  x  =  0.05,  the  strain  response  to  the  applied  load.  The  exact  solution 
is  shown  as  the  thick  black  line,  while  the  FE  solution  is  shown  as  the  thick  grey  line.  The 
difference  between  these  two  solutions — that  is,  the  error  in  the  FE  solution — is  shown  as 
the  thin  black  line. 

As  was  expected,  the  FE  solution  diverges  more  from  the  exact  solution  as  time  pro¬ 
gresses.  The  initial  cycles  (up  to  approximately  1.8  non-dimensional  time  units)  has  a 
comparatively  small  error,  and  hence  we  use  this  portion  of  the  transient  strain  response 
to  predict  the  load  on  the  FE  beam.  Incidently  the  period  of  this  first  cycle  is  approxi¬ 
mately  the  period  of  the  first  mode,  which  demonstrates  the  dominance  of  the  first  mode. 


5.3  Simulation  Results  of  Load  Prediction 

The  results  are  mainly  shown  in  two  forms:  (i)  plots  of  the  predicted  load  coefficients 
versus  mode  number  and  (ii)  comparison  plots  of  the  predicted  load  and  applied  load  along 
the  beam. 

The  load  coefficients  (that  is,  the  ooj  in  Equation  (4.12))  were  determined  using  the 
technique  developed  in  §  4.3  and  the  strains  derived  from  the  transient  FE  solution.  To 
reduce  the  error  introduced  by  the  discretisation  within  the  FE  package  Nastran  (see 
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time  (non-dimensional) 


Figure  5.2:  The  strain  response  at  x  —  0.05  to  a  constant  load  x^(l  —  x) 
applied  as  a  step  function  at  time  t  =  0.  The  exact  and  FE  solutions 
are  shown  as  the  thick  black  and  thick  grey  lines  respectively,  while  the 
error  in  the  FE  solution  is  shown  as  a  thin  black  line. 


Figure  5.2),  the  simulated  strain  measurements  were  taken  exclusively  from  within  the 
first  period  (that  is,  t  ^  1.8). 

Figure  5.3  shows  a  comparison  between  the  predicted  and  the  exact  load  coefficients. 
As  can  be  seen,  the  error  is  larger  in  the  higher  modes  than  the  lower  modes.  The  first  forty- 
nine  simulated  strain  measurements  were  used  to  determine  the  load  coefficients  shown 
in  Figure  5.3.  This  number  of  strain  measurements  was  chosen  because  it  minimised  the 
prediction  error.  We  will  see  later  that  the  number  of  strain  measurements  used  to  predict 
the  load  coefficient  has  a  strong  effect  on  the  resulting  prediction  error. 


Trying  to  predict  the  first  ten  load  coefficients  with  forty-nine  strain  measurement 
resulted  in  an  over-determined  system  of  equations.  We  solved  this  over-determined  system 
using  singular  value  decomposition,  which  is  equivalent  to  using  least  squares  (see,  for 
example,  Golub  and  van  Loan  [10,  p.  242]). 

Figure  5.4  shows  a  plot  of  the  predicted  load  distribution  on  the  FE  beam  model. 
Notice  the  relatively  high  frequency  oscillation  of  the  predicted  load  (broken  grey  line) 
about  the  exact  solution  (solid  black  line).  This  oscillation  is  predicted  by  the  larger 
errors  of  the  higher  modes  shown  in  Figure  5.3. 


In  order  to  measure  how  well  our  load  prediction  technique  works,  we  define  a  nor¬ 
malised  error  in  terms  of  the  coefficients  of  the  load’s  modes 


normalised  error  = 


(5.3) 
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mode  number 

Figure  5.3:  Comparison  of  the  exact  load  coefficients  (denoted  by  Q) 
and  the  predicted  load  coefficients  from  the  FE  derived  strains  ( denoted 
by 

where  N  is  the  number  of  modes  in  our  exact  solution  (for  us  N  =  10),  aoj  are  the 
coefficients  of  the  exact  solution’s  load  modes  (see  Equation  (4.12)),  and  doj  are  our 
predictions  of  oo^.  The  summation  in  the  numerator  goes  up  to  N,  so  if  we  are  only 
predicting  the  first  n  modes  of  the  load  (in  Figure  5.5,  for  example,  n  is  either  three,  six, 
or  nine),  then  ooj  =  0  for  j  =  n  +  1,  n  +  2, . . . ,  A^. 

The  definition  of  normalised  error  given  in  Equation  (5.3)  is  in  fact  closely  related  to 
the  error  in  terms  of  absolute  areas.  This  statement  shouldn’t  come  as  a  surprise  given 
that  the  form  of  this  normalised  error  resembles  Parseval’s  identity  for  the  Fourier  series 
(see,  for  example,  Spiegel  [25]).  In  other  words,  the  normalised  error  is  approximately 
equal  to  the  area  between  the  predicted  and  exact  load  distributions  divided  by  the  area 
under  the  exact  load  distribution  (refer  to  Figure  5.4). 

We  now  justify  the  above  statement.  Substituting  Equation  (4.2)  for  the  load  distri¬ 
bution  we  have  that 


N 

j=i 

where  we  have  made  use  of  Equation  (3.4)  (that  is,  the  orthonormality  of  the  beam  modes) 
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Figure  5-4:  Plot  of  predicted  load  distribution  using  the  transient  strain 
response  of  the  FE  model.  The  load  applied  to  the  FE  model  is  shown 
as  the  solid  black  line,  while  the  predicted  load  is  shown  as  a  broken  grey 
line.  The  horizontal  axis  is  the  non-dimensional  distance  from  the  fixed 
end  of  the  beam. 


in  going  from  the  second  last  line  to  the  last  line.  In  the  first  line  of  the  above  expressions 
the  left  and  right  hand  sides  are  approximately  equal,  instead  of  exactly  equal,  because 
the  number  of  modes  in  the  “exact”  solution  is  finite.  Using  Equation  (5.4)  and  Holder’s 
inequality  for  integrals  (see,  for  example,  Spiegel  [25])  we  have  that 

^1  (  i 

\f{x,t)\dx  <  [f{x,t)fdxj 

The  above  inequality  shows  the  relation  between  the  area  under  the  load  and  the  coeffi¬ 
cients  of  the  load’s  modes. 

As  we’ve  already  noted  the  error  in  predicting  the  load  was  sensitive  to  the  number 
of  strain  measurements  used.  In  addition,  the  load  prediction  error  was  also  sensitive  to 
the  number  of  load  modes  we  were  trying  to  predict.  Figure  5.5  shows  how  error  varies 
with  these  two  parameters.  The  black,  dark  grey,  and  light  grey  lines  represent  the  load 
prediction  of  the  first  three,  the  first  six,  and  the  first  nine  modes  respectively.  The 
horizontal  axis  shows  the  number  of  strain  measurements  used  to  predict  the  load,  while 
the  vertical  axis  shows  the  normalised  error  of  the  resulting  prediction. 

In  order  to  use  Figure  5.5  to  determine  how  well  our  prediction  technique  works,  we 
need  to  know  the  quality  of  the  FE  results  from  which  we  are  deriving  the  measured  strains. 
The  relative  error  in  the  FE  solution’s  strain  response  was  calculated  as:  the  root  mean 
square  (RMS)  of  the  difference  between  the  exact  and  FE  solution  divided  by  the  RMS  of 
the  exact  solution.  For  the  first  fifty  strain  measurements,  the  normalised  error  in  the  FE 
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Figure  5.5:  The  normalised  error  of  the  load  prediction  (vertical  axis)  for 
several  different  configurations.  The  horizontal  axis  shows  the  number 
of  strain  measurements  used  to  predict  the  load.  The  number  of  modes 
used  to  predict  the  load  varied  from  three  o,  through  six  x,  to  nine  k 
modes. 


response  was  approximately  0.03,  that  is,  about  3%  error.  From  Figure  5.5,  we  see  that 
only  the  first  three  modes  could  be  predicted  to  within  about  10%  error,  which  is  larger 
than  the  error  in  the  measurements.  The  error  in  predicting  a  high  number  of  modes, 
nine  for  example,  resulted  in  excessively  large  errors  (as  shown  in  Figure  5.4,  which  was  a 
prediction  example  with  the  least  possible  error).  These  preliminary  results  suggest  that 
the  load  prediction  technique  is  sensitive  to  noise,  and  as  we  show  below,  this  is  indeed 
the  case  for  higher  mode  predictions. 

Before  leaving  Figure  5.5,  note  the  strong  relation  between  the  error  and  the  number 
of  strain  measurements  used  to  predict  the  load.  In  most  cases,  the  optimal  solution  is 
an  order  of  magnitude  more  accurate  than  most  other  solutions  (in  fact,  several  orders  of 
magnitude  better  than  an  ill  chosen  number  of  measurements).  We  also  see  that  the  error, 
after  an  initial  jump,  starts  decreasing  to  its  optimal  value,  after  which  time  the  error  rises 
with  the  number  of  measurements  used  in  the  prediction.  We  can  easily  explain  the  rise 
in  errors  with  increasing  measurements  by  referring  back  to  Figure  5.2,  where  we  see  that 
the  errors  in  the  FE  solution  rise  with  increasing  time. 

Without  knowing  the  exact  solution  a  priori  it’s  hard  to  determine  the  optimal  num¬ 
ber  of  measurements  to  use  in  load  prediction.  As  we  see  from  Figure  5.5  if  we  choose  too 
few  measurements  the  results  could  be  meaningless,  while  choosing  too  many  measure¬ 
ments  leads  us  away  from  the  optimal  solution.  Although  it’s  better  to  choose  too  many 
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measurements  than  too  few,  it’s  best  to  get  it  right.  One  way  to  achieve  a  near  optimal 
solution  would  be  to  plot  the  RMS  of  the  difference  between  successive  predictions  of  load 
coefficients  versus  the  number  of  measurements  used  for  the  prediction.  We  should  see  a 
minimum  near  the  optimal  solution.  In  other  words,  the  coefficients  are  relatively  stable 
near  the  optimal  solution  (see  the  error  of  the  three  mode  solution  in  Figure  5.5).  A  plot 
of  RMS  difference  between  successive  solutions  would  also  alert  the  user  to  unstable  solu¬ 
tions,  which  are  probably  erroneous  (for  example  see  the  errors  of  the  nine  mode  solutions 
in  Figure  5.5). 

We  now  return  to  the  issue  of  load  prediction  sensitivity  to  noise.  In  order  to  determine 
the  effects  of  noise  on  load  prediction,  we  will  add  1%  random  white  noise  to  the  exact 
strain  response.  We  will  then  predict  the  load  distribution  of  this  contaminated  signal. 

Figure  5.6  shows  what  happens  when  we  try  to  predict  the  higher  load  modes  of  the 
exact  strain  response  with  1%  white  noise  added.  Figure  5.7  shows  the  predicted  load 
response.  As  can  be  seen  from  these  plots  it  only  takes  a  small  amount  of  error  to  destroy 
the  higher  mode  signal.  Twenty-five  strain  measurements  were  used  to  derive  the  load 
prediction  results  shown  in  Figure  5.6  and  5.7;  this  number  of  measurements  resulted  in 
the  best  possible  approximation.  Using  either  more  or  fewer  strain  measurements  led  to 
load  prediction  results  with  significantly  larger  errors. 
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Figure  5.6:  Comparison  of  the  exact  load  coefficients  (denoted  by  Q) 
and  the  predicted  load  coefficients  for  the  exact  strain  response  with  1% 
white  noise  added  (denoted  by  13/ 


In  contrast,  the  low  mode  signals  are  quite  robust  as  illustrated  in  Figure  5.8,  which 
shows  the  prediction  of  only  the  lower  load  modes.  Forty-five  strain  measurements  were 
used  to  derive  the  load  prediction  results  shown  in  Figure  5.8,  which  resulted  in  a  load 
prediction  with  the  lowest  error.  Due  to  the  stability  of  the  solution,  choosing  a  larger 
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Figure  5.7:  Plot  of  predicted  load  distribution  along  the  beam  using  the 
transient  exact  strain  response  with  1%  white  noise  added.  The  exact 
applied  load  is  shown  as  the  solid  black  line,  while  the  predicted  load  is 
shown  as  a  broken  grey  line.  The  horizontal  axis  is  the  non-dimensional 
distance  from  the  fixed  end  of  the  beam. 


number  of  strain  measurements  led  to  only  marginally  poorer  load  predictions. 
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Figure  5.8:  Plot  of  predicted  load  distribution  (lower  modes  only)  using 
the  transient  exact  strain  response  with  1%  white  noise  added.  The  exact 
applied  load  is  shown  as  the  solid  black  line,  while  the  predicted  load  is 
shown  as  a  broken  grey  line.  The  horizontal  axis  is  the  non-dimensional 
distance  from  the  fixed  end  of  the  beam. 
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6  Discussion 


Let  us  now  consider  the  assumptions  we  have  made;  we  will  order  these  assumptions 
from  most  restrictive  to  least  restrictive.  Future  directions  for  load  determination  are  then 
outlined. 

The  most  restrictive  assumption  (in  fact  a  requirement)  is  that  we  know  the  initial 
deformation  to  determine  the  loading  distribution  (which  is  constant  in  time).  (We  later 
show  how  to  overcome  the  restriction  of  a  loading  distribution  that  is  constant  in  time.) 
In  reality,  for  continuous  monitoring  of  rotor  blades  we  would  not  have  this  information. 
As  such  this  requirement  of  knowing  the  initial  deformation  places  a  severe  restriction  on 
using  this  approach  to  determine  the  rotor  blade  loading  distribution. 

The  second  most  restrictive  assumption  (which  was  an  implicit  one)  is  the  form  of  the 
PDF  governing  vibration  in  the  beam.  The  implicit  assumption  made  in  using  the  PDE 
given  by  Equation  (2.4)  is  that  the  beam’s  vibration  is  not  damped.  In  reality,  we  would 
expect  the  aerodynamic  damping  to  play  some  role,  and  this  may  cause  coupling  between 
the  modes. 

Likewise  we  have  ignored  chordwise  effects  without  knowing  their  influence  (if  any). 
Quite  clearly  if  the  rotor  blades  twist,  then  the  stiffness  properties  woirld  change  resulting 
in  a  different  vibration  response  (stiffness  variation  is  further  discussed  below). 

The  formulation  we  have  undertaken  has  assumed  that  all  natural  frequencies  (an 
infinite  number  of  them)  are  measured.  In  reality  all  the  infinite  summations  developed 
in  this  report  would  be  truncated  to  finite  summations,  which  would  lead  to  an  aliasing 
effect.  We  stated  earlier  (see  page  5)  that  the  natural  frequencies  of  a  cantilever  beam  are 
given  by  0^^/ETJ{n¥).  Using  Equation  (3.9)  we  then  have  that  the  natural  frequencies 
are 


That  is,  the  natural  frequencies  go  up  quadratically  with  mode  number,  and  hence  it 
becomes  increasingly  harder  to  capture  higher  natural  frequencies.  We  would  expect  that 
the  high  frequency  modes  would  have  small  amplitudes,  and  so  the  aliasing  effects  would 
be  of  second  order  as  compared  to  the  low  frequency  mode  amplitudes.  This  expectation 
was  borne  out  by  the  simulations  in  §  5. 

The  location  of  the  accelerometer  would  also  be  crucial  under  some  circumstances. 
It  would  be  imprudent  to  place  the  accelerometer  at  one  of  the  lower  mode’s  stationary 
nodes.  For  example,  placing  the  accelerometer  at  x  ~  0.358  would  almost  certainly  reduce, 
if  not  completely  eliminate,  the  vibration  signal  from  the  fourth  mode  (see  Figure  3.1). 
Fortunately,  only  the  higher  modes  of  a  cantilever  beam  have  nodes  near  the  fixed  end. 
For  example,  placing  the  accelerometer  at  x  Ri  0.05  (that  is,  5%  along  the  beam  from 
the  fixed  end)  would  only  a  ffect  modes  equal  to  or  higher  than  the  twenty-fifth  mode. 
In  general,  higher  modes  will  have  their  first  stationary  node  at  x  5/(4j  —  2),  where  j 
is  the  mode  number.  This  approximation,  of  the  first  root  of  the  jth  mode,  comes  from 
Equation  (3.6)  and  makes  use  of  the  relation  cosh^  — sinh0  =  e~^  and  the  approximations 
given  by  Equatioirs  (3.9)  and  (3.12). 

Using  the  same  principles,  we  can  extend  this  load  determination  technique  (using 
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vibration  at  a  single  point)  to  non-constant  mass  and  stiffness  problems.  However,  we 
would  then  be  forced  to  solve  the  governing  PDE  numerically,  which  would  add  an  extra 
degree  of  difficulty.  We  have  already  mentioned  that  a  twisting  blade  would  have  varying 
stiffness.  Similarly  centrifugal  effects  would  add  considerable  stiffness  to  a  rotor  blade;  that 
is,  different  rotational  speeds  would  produce  different  blade  stiffnesses.  The  differential 
equation  (DE)  governing  the  vibration  response  of  a  beam  subject  to  an  axial  load  F  is 
given  by  Kelly  [14,  p.  502]  as 


El 


dx'^ 


dx^ 


Notice  the  only  difference  between  the  above  equation  and  the  DE  we  used  (Equation  (2.4)) 
is  the  term  Fd^y/dx^.  The  above  expression  is  derived  by  Meirovitch  [19,  p.  440-3]  from 
first  principles. 

We  have  only  briefly  touched  on  the  effects  of  noise,  especially  the  sensitivity  of  the 
solution  with  respect  to  noise.  Intuitively  we  would  expect  that  noise  would  only  be  a 
problem  for  higher  frequency  modes  (namely  the  modes  we  have  truncated),  and  so  noise 
may  not  be  a  cause  for  concern.  However,  the  results  shown  in  Figure  5.6  (and  associated 
simulations  of  §  5)  cast  doubt  on  this  intuition. 

On  a  helicopter  we  probably  couldn’t  easily  attach  an  accelerometer  or  a  strain  gauge 
to  a  rotor  blade,  and  hence  would  have  to  place  it  on  some  adjacent  structure.  This 
adjacent  structure  may  be  subject  to  damping  and  vibration  itself,  which  may  act  as  noise 
on  the  rotor  blade’s  vibration  signal.  A  similar  but  more  severe  problem  would  arise 
from  the  vibration  of  the  non-measured  rotor  blades.  All  the  blades  would  have  almost 
identical  vibration  responses,  and  hence  it  would  be  hard  (if  not  impossible)  to  distinguish 
the  vibration  signal  from  the  desired  blade  (especially  if  measurements  were  undertaken 
on  some  adjacent  structure). 

In  §  4.3  (the  most  promising  approach)  we  assume  the  loading  could  be  decomposed 
into  the  shape  of  the  natural  modes  (see  Equation  (4.11)).  From  the  boundary  condi¬ 
tions  (Equations  (3.1))  and  Figure  3.1  we  see  that  all  the  natural  modes  have  both  zero 
deformation  and  zero  slope  at  the  fixed  end  of  the  beam  {x  =  0).  Hence  we  cannot 
model  a  non-zero  load  at  the  fixed  end.  In  fact,  non-zero  loads  close  to  the  fixed  end  will 
necessarily  involve  high  frequency  mode  shapes  in  their  representation  (due  to  the  Gibbs 
phenomenon) . 

Superficially,  the  assumption  of  temporally-constant  loading  distribution  appears  very 
severe.  However,  there  is  nothing  to  stop  us  partitioning  the  rotation  of  a  helicopter 
blade  into  several  time  segments,  and  considering  the  loading  constant  over  each  of  these 
segments.  Given  that  the  rotation  frequency  of  the  main  rotor  on  helicopters  is  relatively 
low  (4.3  Hz  for  the  Black  Hawk)  as  compared  to  strain  recording  capabilities,  we  could 
easily  carry  out  this  segmentation  approach  to  reconstruct  a  temporally-varying  loading 
distribution.  Figure  6.1  shows  an  example  of  how  we  can  approximate  a  temporally- 
varying  load  by  discretising  the  load  into  several  steps.  We  see  in  Figure  6.1  that  the  load 
is  constant  over  each  time  interval,  and  changes  in  steps  between  time  intervals. 
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Figure  6.1:  An  example  of  how  a  continuous  time-varying  load  f{t) 
(black  curve)  can  be  discretised  into  steps  (grey  curve),  thus  allowing 
the  use  of  the  assumption  of  a  “constant”  loading  distribution  over  each 
time  interval. 


7  Conclusion 

After  developing  the  non-dimensional  partial  differential  equation  governing  the  vi¬ 
bration  of  a  cantilever  beam,  we  solved  the  homogeneous  form  using  the  separation  of 
variables  technique. 

In  the  main  section  of  this  report  (§  4),  we  showed  that  under  certain  circumstances 
we  could  predict  the  load  distribution  along  a  cantilever  beam.  For  this  load  prediction  we 
needed  the  initial  deformation  of  the  beam  and  either  the  strain  response  or  the  vibration 
response  at  one  location  along  the  beam.  The  reason  we  needed  to  know  the  initial 
deformation  was  that  it  and  the  strain  (or  vibration)  response  were  related,  and  could  not 
be  separated  without  additional  information.  The  need  to  know  the  initial  deformation 
places  a  severe  restriction  on  the  practical  implementation  of  the  load  prediction  technique 
outlined  in  this  report. 

Using  a  finite  element  simulation  we  were  able  to  verify  our  load  prediction  technique. 
As  was  expected,  this  simulation  also  demonstrated  the  effect  noise  might  have  on  load 
prediction.  Namely,  noise  corrupted  the  necessary  information  for  the  higher  load  mode 
determination.  In  contrast,  the  low  load  modes  were  seen  to  be  robust. 

Using  a  single  measurement  location  to  determine  the  load  distribution  along  an  entire 
beam,  knowledge  of  the  initial  beam  deformation  was  shown  to  be  the  most  restrictive  as¬ 
sumption.  In  operating  conditions,  the  initial  deformation  of  a  helicopter  rotor  blade  would 
be  unknown,  making  the  implementation  of  this  load  determination  technique  problematic. 
As  such,  it  is  suggested  that  no  further  consideration  be  given  to  this  load  determination 
approach  using  a  single  measurement  location.  Under  normal  operating  conditions,  multi¬ 
ple  measurement  locations  would  be  required  to  determine  the  loading  distribution  along 
a  rotor  blade. 
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